Set-up

library(here) # for reproducible paths
library(SingleCellExperiment)
library(scater) # For qcs
library(org.Mm.eg.db) # To annotate the genenames
sce <- readRDS(here("processed", "sce_old.RDS"))
age <- "old"

The object has 32285 genes and 23007 cells before filtering

Add cell QC metrics to the sce

First we need to sort the gene names and gene symbols, because the default ensembl notation is not very handy. And then save the mitochondrial genes as such.

if (!file.exists(here("processed", age, "sce_preliminary.RDS"))) {
  # obtain full genenames
  genename <- mapIds(org.Mm.eg.db,
                     keys = rownames(sce),
                     keytype = "ENSEMBL",
                     column = c("GENENAME")
                     )
  # Use the symbols as rownames
  # first make gene names unique
  # TODO: save duplicate gene name list
  symb_unique <- uniquifyFeatureNames(rownames(sce), rowData(sce)[, "Symbol"])
  # Now they can be used as rownames
  rownames(sce) <- symb_unique
  # Add full gene names and the uniuqe symbols to the rowdata
  rowData(sce)$symb_uniq <- symb_unique
  rowData(sce)$gene_name <- genename
  # Subset the mitochondrial genes
  is_mito <- grepl("^mt-", rownames(sce))
} else {
  sce <- readRDS(here( "processed", age, "sce_preliminary.RDS"))
}

Then we can use the scater package to add the quality per cell. This computes for each cell some useful metrics such as the number of umi counts (library size), the number of detected genes and the percentage of mitochondiral genes.

Then we use the automatic isOutlier function from the same package that determine which values in a numeric vector are outliers based on the median absolute deviation (MAD). When using this function with low number a log transformation is added, that prevents negative thresholds.

if (!file.exists(here( "processed", age, "sce_preliminary.RDS"))) {
  sce <- addPerCellQC(sce, subsets = list(mt = is_mito))

  # Automated outlier detection
  outlier_lib_low <- isOutlier(sce$total, log = TRUE, type = "lower")
  outlier_expr_low <-
    isOutlier(sce$detected, log = TRUE, type = "lower")
  outlier_lib_high <- isOutlier(sce$total, type = "higher")
  outlier_expr_high <-
    isOutlier(sce$detected, type = "higher")
  outlier_mt <- isOutlier(sce$subsets_mt_percent, type = "higher")
  # total
  outlier <-
    outlier_lib_low |
      outlier_expr_low |
      outlier_lib_high | outlier_expr_high | outlier_mt

  # Visualize the thresholds and the cells deleted by each parametre
  summary_outlier <- data.frame(
    lib_size_high = c(sum(outlier_lib_high),
                      attr(outlier_lib_high, "thresholds")[2]),
    expression_high = c(sum(outlier_expr_high),
                        attr(outlier_expr_high, "thresholds")[2]),
    lib_size_low = c(sum(outlier_lib_low),
                     attr(outlier_lib_low, "thresholds")[1]),
    expression_low = c(sum(outlier_expr_low),
                       attr(outlier_expr_low, "thresholds")[1]),
    mt_pct = c(sum(outlier_mt), attr(outlier_mt, "thresholds")[2]),
    total = c(sum(outlier), NA)
  )
  row.names(summary_outlier) <- c("Cells filtered", "Threshold")
  write.csv(summary_outlier, here("outs", age, "autofilter_summary.csv"))
  # Add if it is an outlier to the metadata
  sce$outlier <- outlier
} else {
  summary_outlier <- read.csv(here("outs", age, "autofilter_summary.csv"))
}
summary_outlier
##                X lib_size_high expression_high lib_size_low expression_low
## 1 Cells filtered        956.00          93.000       0.0000       773.0000
## 2      Threshold      16021.98        5867.663     314.1608       162.4078
##       mt_pct total
## 1 2724.00000  4381
## 2   21.74347    NA

Plots before QC

Diagnostic plots to visualize the data distribution. The orange cells are marked as outliers by the automatic detection from scater. These thresholds will probably be manually adjusted to better fit the distribution of our data.

Violin plots

plotColData(sce, x = "Sample", y = "sum", colour_by = "outlier") +
  ggtitle("Total count") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

plotColData(sce, x = "Sample", y = "sum", colour_by = "outlier") +
  scale_y_log10() + ggtitle("Total count log scale") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

plotColData(sce, x = "Sample", y = "detected", colour_by = "outlier") +
  scale_y_log10() + ggtitle("Detected Genes") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

plotColData(sce, x = "Sample", y = "sum", colour_by = "chip") +
  scale_y_log10() + ggtitle("total count by batch") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

plotColData(sce, x = "Sample", y = "subsets_mt_percent", colour_by = "outlier") +
  ggtitle("Mitocchondrial percentatge") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Histograms

In the x axis we can see the total number of umi (library size) per cell, the number of detected genes per cell and the mitochondrial percentage per cell; with the number of cells for each measure in the y axis.

hist(
  sce$total,
  breaks = 100
)

This object had already been filtrated with the cell-calling algorithm from CellRanger, that is meant to remove empty droplets. Therefore it is expected to see the total sum of umi skewed as in the plot above.

hist(
  sce$detected,
  breaks = 100
)

A bimodal plot is not something usual for the detected number of genes. This bimodality was already visible in the violin plots.

hist(
  sce$subsets_mt_percent,
  breaks = 100
)

There is a very heavy tail of cells with high mitochondrial genes.

Scatter plots

plotColData(sce, x = "sum", y = "subsets_mt_percent", colour_by = "outlier")

plotColData(sce, x = "sum", y = "detected", colour_by = "outlier")

plotColData(sce, x = "sum", y = "detected", colour_by = "Sample")

PCA

Here we run a PCA using the information in the metadata instead of the gene expression. It is useful to visualize the QC parametres.

if (!file.exists(here( "processed", age, "sce_preliminary.RDS"))) {
  sce <- runColDataPCA(sce, variables = c("sum", "detected", "subsets_mt_percent"))
}
plotReducedDim(sce, dimred = "PCA_coldata", colour_by = "Sample")

sce$chip <- as.character(sce$chip)
plotReducedDim(sce, dimred = "PCA_coldata", colour_by = "chip")

Ratio between sum and gene counts

This measures the number of detected genes per cell dividided by its lirary size. This will be very useful to delete the cells that have low gene counts but a relatively high umi count.

if (!file.exists(here( "processed", age, "sce_preliminary.RDS"))) {
  sce$ratio_detected_sum <- sce$detected / sce$sum
  sce$outlier_ratio <- isOutlier(sce$ratio_detected_sum, type = "both")
}

summary(sce$ratio_detected_sum)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.005707 0.339102 0.393226 0.394572 0.455339 0.828814
plotColData(sce, x = "sum", y = "detected", colour_by = "ratio_detected_sum")

# Plot an histogram with the ratio between umi and gene counts
hist(
  sce$ratio_detected_sum,
  breaks = 100
)

The isoutlier function can but used to find the outliers of any distribution, as far as it is roughly normal. Bellow we use it with the ratio between the number of genes expressed and the number of umi.

# Use the is outlier function from scater to see the cutoffs suggestions
plotColData(sce, x = "sum", y = "detected", colour_by = "outlier_ratio")

attr(sce$outlier_ratio, "thresholds")
##     lower    higher 
## 0.1365831 0.6498693

The suggestions here are very reasonable, it is not enough to delete all the cells we do not trust but we can definitely take this into consideration as a threshold parameter.

Decide the thresholds

The upper thresholds from the sum of umi counts, that allows filtering doublet cells, might be a bit too harsh. This could be due to the fact the distribution was not normal. We could start with 25,000 instead of 1.6022^{4}.

The other thresholds seem more reasonable. These are deleting cells with less than 5868 or more than 162 detected genes, and less than 314 umi counts.

Finally, we will also take in consideration the detected genes/umi counts ratio, that filter out cells with relatively high umi counts but very few detected genes ( lots of copies from the same genes?). As well as the automatic mitochondrial threshold, 21.74 %, this is a bit high compared with other datasets but it is already deleting 2724 cells. Lower down to only 10 % would delete 7297 cells.

Here we decide the new thersholds.

if (!file.exists(here( "processed", age, "sce_preliminary.RDS"))) {
  # New thresholds
  discard_lib_high <- sce$sum > 25000
  discard_expr_high <- outlier_expr_high
  discard_expr_low <- outlier_expr_low
  discard_lib_low <- outlier_lib_low
  discard_ratio <- sce$outlier_ratio
  discard_mt <- outlier_mt
  sce$discard <- discard_lib_high | discard_expr_high | discard_expr_low |
    discard_lib_low | discard_ratio | discard_mt

  # Create a summary with the number of cells deleted and the thresholds

  summary_discard <- data.frame(
    lib_size_high = c(sum(discard_lib_high), 25000),
    expression_high = c(sum(discard_expr_high),
                        attr(outlier_expr_high, "thresholds")[2]),
    lib_size_low = c(sum(discard_lib_low),
                     attr(outlier_lib_low, "thresholds")[1]),
    expression_low = c(sum(discard_expr_low),
                       attr(outlier_expr_low, "thresholds")[1]),
    mt_pct = c(sum(discard_mt), attr(outlier_mt, "thresholds")[2]),
    total = c(sum(sce$discard), NA)
  )

  row.names(summary_discard) <- c("Cells filtered", "Threshold")
  write.csv(summary_discard, here("outs", age, "reason_for_discard_QC.csv"))
  # Save the object at this stage with the label "preliminary analysis"
  # As only annotation addition has been done, without deleting anything yet.
  saveRDS(sce, here( "processed", age, "sce_preliminary.RDS"))
} else {

}
## NULL
summary_discard <- read.csv(here("outs", age, "reason_for_discard_QC.csv"))
summary_discard
##                X lib_size_high expression_high lib_size_low expression_low
## 1 Cells filtered           196          93.000       0.0000       773.0000
## 2      Threshold         25000        5867.663     314.1608       162.4078
##       mt_pct total
## 1 2724.00000  4383
## 2   21.74347    NA

QC with new thresholds

Bellow we are going to plot the same plots as before with the new thresholds

Violin plots

plotColData(sce, x = "Sample", y = "sum", colour_by = "discard") +
  ggtitle("Total count") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

plotColData(sce, x = "Sample", y = "sum", colour_by = "discard") +
  scale_y_log10() + ggtitle("Total count log scale") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

plotColData(sce, x = "Sample", y = "detected", colour_by = "discard") +
  scale_y_log10() + ggtitle("Detected Genes") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

plotColData(sce, x = "Sample", y = "subsets_mt_percent", colour_by = "discard") +
  ggtitle("Mitocchondrial percentatge") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Scatter plots

plotColData(sce, x = "sum", y = "subsets_mt_percent", colour_by = "discard")

plotColData(sce, x = "sum", y = "detected", colour_by = "discard")

Gene QC

It is typically a good idea to remove genes whose expression level is considered “undetectable”. Here we define a gene as detectable if at least two cells contain a transcript from the gene. It is important to keep in mind that genes must be filtered after cell filtering since some genes may only be detected in poor quality cells.

genes_beforeqc <- dim(sce)[1]
keep_feature <- rowSums(counts(sce) > 0) > 1
sce <- sce[keep_feature, ]
genes_beforeqc - dim(sce)[1]
## [1] 9604

This way we deleted 9604 genes and kept 22681 genes

We can look at a plot that shows the top 50 (by default) most-expressed features. Each row in the plot below corresponds to a gene; each bar corresponds to the expression of a gene in a single cell; and the circle indicates the median expression of each gene, with which genes are sorted. We expect to see the “usual suspects”, i.e., mitochondrial genes, actin, ribosomal protein, MALAT1. A large number of pseudo-genes or predicted genes may indicate problems with alignment.

plotHighestExprs(sce, exprs_values = "counts")

Create new filetered object

if (!file.exists(here( "processed", age, "sce_QC.RDS"))) {
  sce <- sce[, sce$discard == FALSE]
  saveRDS(sce, here( "processed", age, "sce_QC.RDS"))
}

The object has 22681 genes and 18624 cells after filtering

Session Info

Click to expand 
sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] org.Mm.eg.db_3.12.0         AnnotationDbi_1.52.0       
##  [3] scater_1.18.6               ggplot2_3.3.3              
##  [5] SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
##  [7] Biobase_2.50.0              GenomicRanges_1.42.0       
##  [9] GenomeInfoDb_1.26.2         IRanges_2.24.1             
## [11] S4Vectors_0.28.1            BiocGenerics_0.36.0        
## [13] MatrixGenerics_1.2.1        matrixStats_0.58.0         
## [15] here_1.0.1                 
## 
## loaded via a namespace (and not attached):
##  [1] viridis_0.5.1             sass_0.3.1               
##  [3] BiocSingular_1.6.0        bit64_4.0.5              
##  [5] jsonlite_1.7.2            viridisLite_0.3.0        
##  [7] DelayedMatrixStats_1.12.3 scuttle_1.0.4            
##  [9] bslib_0.2.4               highr_0.8                
## [11] blob_1.2.1                GenomeInfoDbData_1.2.4   
## [13] vipor_0.4.5               yaml_2.2.1               
## [15] pillar_1.5.1              RSQLite_2.2.3            
## [17] lattice_0.20-41           glue_1.4.2               
## [19] beachmat_2.6.4            digest_0.6.27            
## [21] XVector_0.30.0            colorspace_2.0-0         
## [23] htmltools_0.5.1.1         Matrix_1.3-2             
## [25] pkgconfig_2.0.3           zlibbioc_1.36.0          
## [27] scales_1.1.1              BiocParallel_1.24.1      
## [29] tibble_3.1.0              farver_2.1.0             
## [31] ellipsis_0.3.1            cachem_1.0.4             
## [33] withr_2.4.1               magrittr_2.0.1           
## [35] crayon_1.4.1              memoise_2.0.0            
## [37] evaluate_0.14             fansi_0.4.2              
## [39] beeswarm_0.3.1            tools_4.0.4              
## [41] lifecycle_1.0.0           stringr_1.4.0            
## [43] munsell_0.5.0             DelayedArray_0.16.2      
## [45] irlba_2.3.3               compiler_4.0.4           
## [47] jquerylib_0.1.3           rsvd_1.0.3               
## [49] rlang_0.4.10              grid_4.0.4               
## [51] RCurl_1.98-1.2            BiocNeighbors_1.8.2      
## [53] bitops_1.0-6              labeling_0.4.2           
## [55] rmarkdown_2.7             gtable_0.3.0             
## [57] DBI_1.1.1                 R6_2.5.0                 
## [59] gridExtra_2.3             knitr_1.31               
## [61] fastmap_1.1.0             bit_4.0.4                
## [63] utf8_1.1.4                rprojroot_2.0.2          
## [65] stringi_1.5.3             ggbeeswarm_0.6.0         
## [67] Rcpp_1.0.6                vctrs_0.3.6              
## [69] xfun_0.21                 sparseMatrixStats_1.2.1